{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xJWpFX-zc-fE"
   },
   "source": [
    "# Hydrostats Features\n",
    "Hydrostats is a module for time series comparison. It has preproccessing tools, visualization tools, a library of over 50 error metrics for quantitative analysis of data from two timeseries, and analysis tools. In this notebook examples are given for most of the hydrostats features. The two code blocks below shows how to install hydrostats with PIP and show the current version (make sure it is the most recent version). It also imports the necessary modules for this demonstration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Y6-0WTL9fLMt"
   },
   "outputs": [],
   "source": [
    "# Run these commands if you have not yet installed hydrostats, and make sure it is the latest version\n",
    "\n",
    "# !pip install hydrostats --upgrade\n",
    "!pip show hydrostats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "bFbDkj76NfLn"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "pd.options.display.max_rows = 10\n",
    "import hydrostats as hs\n",
    "import hydrostats.data as hd\n",
    "import hydrostats.visual as hv\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "U_y9uoVydvPD"
   },
   "source": [
    "## Preprocessing\n",
    "Hydrostats has tools that simply the preprocessing of data so that an analysis can be run on two timeseries. Once the data has been preprocessed then the analysis can be run on the two time-series. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "LxXZx4IQeXNF"
   },
   "source": [
    "### Merge Data Function\n",
    "The first thing we need to do is load our data into the notebook, which is done below. In this example we will be comparing forecasts from the Streamflow Prediction Tool (SFPT) Tethys app to forecasts from the Global Flood Awareness System (GLOFAS). The data is loaded from the Hydrostats github repository, but if a user wanted to load local files, they would need to supply the path to the files instead. In this example the data is daily for both datasets. The merge_data Hydrostats function function used below allows a user to take two timeseries and merge them into a combined pandas dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 424
    },
    "colab_type": "code",
    "id": "vrpOHXBhjM4q",
    "outputId": "0e173ed6-4da9-4d81-bc46-d1afa2a271bc"
   },
   "outputs": [],
   "source": [
    "# Defining the URLs of the datasets\n",
    "github_base_url = r'https://raw.github.com/BYU-Hydroinformatics/Hydrostats/master/Sample_data'\n",
    "sfpt_url = github_base_url + r'/sfpt_data/magdalena-calamar_interim_data.csv'\n",
    "glofas_url = github_base_url + r'/GLOFAS_Data/magdalena-calamar_ECMWF_data.csv'\n",
    "\n",
    "merged_df = hd.merge_data(sfpt_url, glofas_url, column_names=['SFPT', 'GLOFAS'])\n",
    "merged_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "SPgxO_us6C6z"
   },
   "source": [
    "### Seasonal Periods\n",
    "Once we have the entire dataset merged into one combined dataframe, can find then go on to find certain seasonal periods that we want to analyze with the seasonal_period hydrostats function. In the example below, we only want to analyze from April 1st to July 31st seasonally within the time range of Jan. 1st, 1986 to Dec. 31st, 1992."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 424
    },
    "colab_type": "code",
    "id": "y9x3hzQSfE07",
    "outputId": "70092e22-4746-44bc-e5a8-ab181cf1c21e"
   },
   "outputs": [],
   "source": [
    "seasonal_df = hd.seasonal_period(merged_df, ['04-01', '07-31'], time_range=['1986-01-01', '1992-12-31'])\n",
    "seasonal_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "tG2ctFQKIKVU"
   },
   "source": [
    "### Daily and Monthly Averages\n",
    "We can also easily find the daily and monthly averages in streamflows over the 35 years of daily data we have. In the example below we use the daily_average function, but a monthly_average function() in the same module is also available. We create a new dataframe that contains the daily averages for this time-series below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 424
    },
    "colab_type": "code",
    "id": "UelqZSJ1Ikcr",
    "outputId": "7eebd3bf-f258-4e3b-e79a-8923703cb195"
   },
   "outputs": [],
   "source": [
    "daily_avg_df = hd.daily_average(merged_data=merged_df)\n",
    "daily_avg_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Pe_tNa4AIuy9"
   },
   "source": [
    "### Daily and Monthly Standard Deviation and Standard Error\n",
    "Hydrostats also can calculate daily and monthly standard deviation as well as daily and monthly standard error. We create a new dataframe below that contains the daily standard error for this time-series below.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 424
    },
    "colab_type": "code",
    "id": "8NIS5O4TJXop",
    "outputId": "59ed2eb7-968c-4df6-ca76-78c1353162e5"
   },
   "outputs": [],
   "source": [
    "daily_std_error = hd.daily_std_error(merged_data=merged_df)\n",
    "daily_std_error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ynITzPTYhdZ5"
   },
   "source": [
    "## Visualization\n",
    "Once preprocessing on the data is finished, the user can then start to visualize the data with the Hydrostats visual module. There are 4 functions included in the Hydrostats visual model, which are: plot, hist, scatter, and qqplot. These functions are demonstrated below. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FMcv7Y5cihv7"
   },
   "source": [
    "### The Plot Function\n",
    "The plot function will simply take the observed and simulated data and compare them by plotting them on the vertical axis with respect to time. When comparing stream flow rates, this is known as a hydrograph, but the plot function can be used for any time series if desired by the user. An few examples of using the plot function is shown below. By using this function, we are able to plot the entire time series, the seasonal period that we specified, the daily averages, and the gaily averages with error bars. \n",
    "\n",
    "**Note**: If you are viewing these examples in Google's Colaboratory, they will not display properly. You will need to use jupyter notebooks or jupyter lab to view them properly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Cu2caqac8byW"
   },
   "source": [
    "#### Plotting the Entire Time Series\n",
    "The entire time series is plotted below with the hydrostats visual plot function. In this example, the linestyles are changed from the default, a title, legend, and labels are added, and metrics are calculated and displayed on the side. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1588
    },
    "colab_type": "code",
    "id": "UXgmzFE_Jesi",
    "outputId": "ab2cbda0-e8b2-48af-e1a3-de64d5b54af0"
   },
   "outputs": [],
   "source": [
    "import hydrostats.visual as hv\n",
    "\n",
    "hv.plot(merged_data_df=merged_df, \n",
    "        title='Hydrograph of Entire Time Series',\n",
    "        linestyles=['r-', 'k-'], \n",
    "        legend=['SFPT', 'GLOFAS'], \n",
    "        labels=['Datetime', 'Streamflow (cfs)'], \n",
    "        metrics=['ME', 'NSE', 'SA'],\n",
    "        fig_size=(14, 8),\n",
    "        text_adjust=(-0.25, 0.8)\n",
    "       )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "SnfZNt3LNfML"
   },
   "source": [
    "#### Plotting Seasonal Averages with Error Bars\n",
    "Below, it is demonstrated how to plot daily averages with error bars. You simply need to pass in the daily average dataframe that we created earlier in the pre-processing step (daily_avg_df()) and pass in the daily standard error dataframe in the ebars argument (daily_std_error()). We change the linestyle and the color of the ebars below and a title, legend, and labels are added. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1494
    },
    "colab_type": "code",
    "id": "9SGghzNfGGNE",
    "outputId": "fa46763a-6b88-43fb-feac-0af2281d7436"
   },
   "outputs": [],
   "source": [
    "hv.plot(merged_data_df=daily_avg_df,\n",
    "        title='Daily Average Streamflow (Standard Error)',\n",
    "        legend=['SFPT', 'GLOFAS'],\n",
    "        x_season=True,\n",
    "        labels=['Datetime', 'Streamflow (csm)'],\n",
    "        linestyles=['r-', 'k-'],\n",
    "        fig_size=(14, 8),\n",
    "        ebars=daily_std_error,\n",
    "        ecolor=['r', 'k']\n",
    "        )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "kMQGdRCPNfMQ"
   },
   "source": [
    "### Plotting a Histogram\n",
    "To plot a histogram, we can either use the merged dataframe that we created or two arrays of simulated and observed data. In this case, we use the merged dataframe because we already have it saved to memory. We specify the number of bins and then add a legend, title, and labels. We can compare the distribution of the data with this tool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1312
    },
    "colab_type": "code",
    "id": "ifiL8kyJNfMR",
    "outputId": "6ca44d57-4815-4dd5-affe-9505e542286a"
   },
   "outputs": [],
   "source": [
    "hv.hist(merged_data_df=merged_df,\n",
    "        num_bins=100,\n",
    "        title='Histogram of Streamflows',\n",
    "        legend=['SFPT', 'GLOFAS'],\n",
    "        labels=['Bins', 'Frequency'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "kIljQrKGNfMW"
   },
   "source": [
    "### Plotting a Quantile-Quantile (qq) Plot\n",
    "Similar to the histogram, the qq plot is a means to check whether two datasets come from the same distribution. In the example below, we use the merged dataframe (merged_df) to supply the data. We can also choose to use two array of both simulated and observed data if desired. We then create a title, legend, and labels. We also change the size of the figure (inches)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 857
    },
    "colab_type": "code",
    "id": "ctKW0p4wNfMX",
    "outputId": "a638b2f5-1cde-4910-d44e-e9f1eb0f6b38"
   },
   "outputs": [],
   "source": [
    "hv.qqplot(merged_data_df=merged_df, title='Quantile-Quantile Plot of Data',\n",
    "          xlabel='SFPT Data Quantiles', ylabel='GLOFAS Data Quantiles', legend=True, \n",
    "          figsize=(8, 6))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "IqvMttWDNfMc"
   },
   "source": [
    "### Plotting a Scatter Plot\n",
    "We can plot scatter plots using the scatter() function. In the example below we pass in the data in the form of arrays to demonstrate how to do this. We plot the data in a scatter plot with both a normal scale and a log-log scale. Typically, the Log-Log scale spreads out the data for a more clear representation. In this example, the data looks about the same with either scale. We add a title, labels, and a legend for the lines added.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1553
    },
    "colab_type": "code",
    "id": "Ae76xjX-NfMd",
    "outputId": "f3f0061b-ed18-4157-df44-aaa722b5c74d"
   },
   "outputs": [],
   "source": [
    "SFPT_array = merged_df.iloc[:, 0].values\n",
    "GLOFAS_array = merged_df.iloc[:, 1].values\n",
    "\n",
    "hv.scatter(sim_array=SFPT_array, obs_array=GLOFAS_array, grid=True, title='Scatter Plot (Normal Scale)', \n",
    "           labels=['SFPT', 'GLOFAS'], best_fit=True)\n",
    "plt.show()\n",
    "\n",
    "hv.scatter(sim_array=SFPT_array, obs_array=GLOFAS_array, grid=True, title='Scatter Plot (Log-Log Scale)', \n",
    "           labels=['SFPT', 'GLOFAS'], line45=True, log_scale=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "PI_HxSjbNfMj"
   },
   "source": [
    "## Error Metric Calculations\n",
    "AFter we have done the preprocessing and the visualization, we can do a more in depth analysis of the data by creating tables of metrics and doing a time lag analysis to check for any potential timing errors in the data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "B2wE_9AZ6-GI"
   },
   "source": [
    "### Calculating Error Metric with the HydroErr package\n",
    "\n",
    "Hydrostats has a module that pulls from the HydroErr package that allows users to access over 60 different error metrics. These metrics have data checks to make sure that the data that you are inputing to the metrics does not have NaN or Inf values, and can optionally check for zero and negative values. Below we demonstrate how to calculate error metrics, and highlight how to remove zero values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "eSRgtZng67oB",
    "outputId": "f2978bce-bde5-43b6-a469-f5629f0fa655"
   },
   "outputs": [],
   "source": [
    "# Putting the data from the dataframe into 1D arrays\n",
    "sfpt_array = merged_df.iloc[:, 0].values\n",
    "glofas_array = merged_df.iloc[:, 1].values\n",
    "\n",
    "# Calculating the KGE metric\n",
    "hs.kge_2012(sfpt_array, glofas_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 121
    },
    "colab_type": "code",
    "id": "_OXgr2r78Joh",
    "outputId": "d61bf98d-8c43-4e6e-910f-837806742f0a"
   },
   "outputs": [],
   "source": [
    "# Note that when we remove zeros the indices where zeros are located are given through a warning\n",
    "hs.kge_2012(sfpt_array, glofas_array, remove_zero=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "6LhxSDy2NfMj"
   },
   "source": [
    "### Creating a Table of Error Metrics\n",
    "\n",
    "To create a table of metrics, we use the make_table() function, located in the main hydrostats import. In this function we can specify the time periods we want to analyze (the entire time series will be analyzed by default) as well as the metrics that we want to use. Below we specify the merged dataframe, the metrics we want to use, the time periods to want to evaluate (high and low periods in this case), to remove zeros and negatives, and the location of the station to put in the table. This may raise a lot of errors if the two timeseries have values that need to be removed, so if you would like to ignore the errors, the warnings.simplefilter('ignore') can be applied. \n",
    "\n",
    "To find what abbreviations correspond to which metric, a table is provided in the docs [here](https://waderoberts123.github.io/Hydrostats/ref_table.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1352
    },
    "colab_type": "code",
    "id": "AJ7tZtFRNfMk",
    "outputId": "9e8ff0da-fe31-476b-f910-f569b5b4d811"
   },
   "outputs": [],
   "source": [
    "hs.make_table(merged_dataframe=merged_df, \n",
    "              metrics=['MAE', 'r2', 'ACC', 'NSE', 'KGE (2012)', 'SA'], \n",
    "              seasonal_periods=[['01-01', '03-31'], ['04-01', '06-30'], ['07-01', '09-30'], ['10-01', '12-31']], \n",
    "              remove_neg=True, remove_zero=True, \n",
    "              location='Magdalena')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ensemble Adjusted Error Metrics and Skill Score\n",
    "\n",
    "Hydrostats also includes many common forecast ensemble metrics such as the ensemble mean squared error (EMSE), the continuous ranked probability score (CRPS), the area under the receiver operating characteristic curve (AUROC), and the ensemble adjusted brier score (EABS). Users can also calculate skill scores for all of these metrics comparing their model’s forecast performance to that of a benchmark forecast.\n",
    "\n",
    "Below we give an example of calculating the crps scores for an ensemble forecast and a benchmark and then computing the skill score to compare the two forecasts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hydrostats.ens_metrics as em \n",
    "\n",
    "# Generate random data as an example \n",
    "np.random.seed(6543934)\n",
    "ens_array = (np.random.rand(15, 52) + 1) * 100  # Matrix 15 x 52, range [0, 100)\n",
    "obs_array = (np.random.rand(15) + 1) * 100 # 1D matrix of length 15, range [0, 100)\n",
    "\n",
    "# Compute crps metric values  \n",
    "crps_forecast = em.ens_crps(obs_array, ens_array)['crps'] # Only store the scores in our variable\n",
    "\n",
    "# Compute example benchmark forecast scores\n",
    "crps_benchmark = em.ens_crps(obs_array, ens_array[:, 1:3])['crps'] # Only use first two ensembles as a benchmark\n",
    "\n",
    "# Compute the skill score\n",
    "skill_score = em.skill_score(scores=crps_forecast, bench_scores=crps_benchmark, perf_score=0)\n",
    "print(skill_score)\n",
    "\n",
    "# In this case, using the 52 ensembles is slightly more skillful than using just the first two. \n",
    "# This is random though so this does not mean much."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Hydrostats_Features.ipynb",
   "provenance": [],
   "toc_visible": true,
   "version": "0.3.2"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
